/*
    Copyright 2006-2011 Patrik Jonsson, sunrise@familjenjonsson.org

    This file is part of Sunrise.

    Sunrise is free software; you can redistribute it and/or modify
    it under the terms of the GNU General Public License as published by
    the Free Software Foundation; either version 3 of the License, or
    (at your option) any later version.

    Sunrise is distributed in the hope that it will be useful,
    but WITHOUT ANY WARRANTY; without even the implied warranty of
    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
    GNU General Public License for more details.

    You should have received a copy of the GNU General Public License
    along with Sunrise.  If not, see <http://www.gnu.org/licenses/>.

*/

/** \file

    Declaration of the xfer, ray_stack_object and propagation_history
    classes. */

// $Id$

#ifndef __xfer__
#define __xfer__

#include "config.h"
#include "mcrx-types.h"
#include "ray.h"
#include "optical.h"
#include <vector>
#include <algorithm> 
#include <memory>
#include "hpm.h"
#include <fstream> 
#include "mcrx-debug.h"
#include "boost/lambda/lambda.hpp"
#include "boost/lambda/bind.hpp"
#include "constants.h"
#include "threadlocal.h"
#include "emission.h"

namespace mcrx {
  template<typename, typename> class ray_stack_object;
  template<typename> class propagation_history;
  template <typename, typename, typename> class xfer;
  // and some forward declarations...
  template <typename> class emergence;
}


/** Defines the objects on the ray stack, with the necessary
    information to continue tracing a ray.*/
template <typename T_ray, typename T_cell_tracker> 
class mcrx::ray_stack_object { 
public: 
  /// Indicates whether the ray was pushed before or after scattering.
  int resume_with_; 
  T_cell_tracker cell_;  ///< The cell the ray is located in.
  T_ray ray_; ///< The ray data.

  ray_stack_object (int rw, const T_ray& r, T_cell_tracker c) :
    resume_with_(rw), cell_(c), ray_(r) {};
};

/** Stores history of propagation through cells when using forced
    scattering.  The propagation_history class contains a path length,
    an accumulated optical depth and a pointer to the cell traversed,
    which is put into a vector when the ray is propagated through the
    cells.  In this way we can quickly find where the forced
    scattering occurs. */
template <typename T_cell_tracker>
class mcrx::propagation_history {
private:
  T_densities n_; ///< Traversed column density.
  T_float length_; ///< Accumulated path length.
  T_cell_tracker cell_; ///< Pointer to the grid cell traversed.
public:
  propagation_history (const T_densities& n, T_float l, const T_cell_tracker& c):
    n_(n), length_(l), cell_ (c){};
  const T_densities& n () const {return n_;};
  T_float length () const {return length_;};
  const T_cell_tracker& cell () const {return cell_;};
  /// Comparison operator so that we can use binary search.
  bool operator< (const mcrx::propagation_history<T_cell_tracker>& rhs) const {
    return length() < rhs.length();};
};

/** The xfer class is responsible for performing the ray tracing.  It
    relies on the emission object to create rays, uses the grid object
    to step the ray grid cell by grid cell, and uses the emergence
    object to register the rays leaving the volume. This class can run
    in both monochromatic and polychromatic mode, depending on the
    dust_model_type set. The details of the algorithm used to trace
    the rays is described in the Sunrise papers. */
template <typename dust_model_type, typename grid_type, 
	  typename T_rng_policy = mcrx::global_random>
class mcrx::xfer {
public:
  typedef dust_model_type T_dust_model;
  typedef grid_type T_grid;
  typedef mcrx::T_float T_float;
  typedef typename T_dust_model::T_scatterer T_scatterer;
  typedef typename T_dust_model::T_lambda T_lambda;
  typedef typename T_dust_model::T_biaser T_biaser;
  typedef ray<T_lambda> T_ray;
  typedef emission<typename T_dust_model::T_chromatic_policy::T_base_policy, 
			    T_rng_policy> T_emission;
  typedef emergence<T_lambda> T_emergence;
  typedef typename T_grid::T_cell T_cell;
  typedef typename T_grid::T_cell::cell_tracker T_cell_tracker;
  typedef propagation_history<T_cell_tracker>  T_propagation_history;

private:

  /// The Grid
  T_grid& g;

  /// Pointer to the emergence object, collecting emerging rays. 
  T_emergence& emergence_;
  /// Local emergence object, if keep_local_emergence is true.
  std::auto_ptr<T_emergence> local_emergence_;
  /** Pointer to either the global or the local emergence, whichever
      is used, so we don't have to worry about which one it is in the
      algorithm. */
  T_emergence* current_emergence_;

  /** The emission object emits rays.  Note that we make a COPY of the
      object here to ensure each thread has its own. The emission
      classes have to be reasonable about making shallow copies to
      large chunks of data. \todo Do we really need a copy here? There
      should be no problem making the emission class reentrant. */
  std::auto_ptr<T_emission> emission_;

  /** The dust_model object holds the scatterers and is responsible
      for translating the densities in the grid to optical depth. */
  const T_dust_model& dust_model;

  /// \name Ray tracing variables
  ///@{
  T_cell_tracker current_cell; ///< The grid_cell the ray is currently in.
  T_float scattering_reftau_; ///< Optical depth where scattering will happen.
  const T_float I_last_scatter; ///< Smallest ray intensity we care about.
  T_ray ray_; ///< The ray currently being traced.
  T_biaser b; ///< The biaser object used for the current ray.
  ///@}

  /** \name Cache variables
      To avoid repeatedly allocating memory for temporary arrays in
      shoot and propagate, we make them class members. (And actually
      they are made to point to the same 2D array so they are also
      localized in memory.)
  */
  ///@{
  /// The intensity normalization of the current ray.
  T_lambda norm_;
  /// Incoming intensity in the grid cell when adding intensities in shoot().
  T_lambda i_in_;
  /// Outgoing intensity in the grid cell when adding intensities in shoot().
  T_lambda i_out_;
  /// Optical depth through the grid cell when adding intensities in shoot().
  T_lambda dtau_;
  /// The optical depth to grid edge when forced scattering.
  T_lambda tau_exit_;
  /// (1-exp(-tau_exit)) for forced scattering.
  T_lambda onemexptau_exit_;
  /// The optical depth to the point of scattering.
  T_lambda scattering_tau_;
  /// The original intensity of the ray, saved for restoring later.
  T_lambda original_intensity_;
  /// The intensity to the camera in emerge_from_point.
  T_lambda intensity_to_camera_;
  /// Column density calculated in propagate().
  T_densities dn_; 
  ///@}

  /** \name Forced scattering variables.
    These keep track of the history of the ray.
  */
  ///@{
  /** The propagation_history vector, keeping track of which cells
      have been traversed since last event. */
  std::vector<T_propagation_history > traversed;
  const T_float I_forced; ///< Smallest ray intensity to force scattering. 
  const T_float N_forced; ///< Number of scatterings to be forced.
  ///@}

  const T_float ray_max_i_;
  const int n_scat_min_; ///< Min number of scatterings to register
  const int  n_scat_max_; ///< Max number of scatterings to register
  /** True if intensity information should be registered in the
      grid. */
  bool add_intensities_;
  /** The ray stack object holds rays that have been split and need to
      be processed. */
  std::vector<ray_stack_object<T_ray, T_cell_tracker> > ray_stack_;

  /// Random number generator policy is used for getting random numbers. 
  T_rng_policy rng_;
  
  void emit();// Emits a ray.
  // Propagates the ray to the boundary of the grid cell. 
  template <bool do_scatter, bool keep_traversed, bool add_intensity>
  bool propagate(T_float max_length= blitz::huge(T_float()));

  void scatter(int resume_with = -1);   // Scatters ray.
  // Calculates emerging intensity from current ray position.
  void emerge_from_point (const angular_distribution<T_lambda>&, const vec3d&);
  // Propagates an external ray to the grid boundary.
  void propagate_external (T_float max_length = blitz::huge(T_float())); 

  /** Adds the ray to the image. */
  template <typename T>
  void add_to_camera(typename T_emergence::T_camera* c, const vec3d& pos, 
		     const T& intensity, T_float doppler) {
    if ((ray_.scatterings()<=n_scat_max_)&&(ray_.scatterings()>=n_scat_min_))
      c->add(pos, intensity, doppler);
  };

  void add_intensity (const T_cell_tracker&, T_float);
  template <typename T>
  void add_intensity (const T_cell_tracker&, T_float, const blitz::ETBase<T>& intensity);

  void check_ray_split(int resume_with=-1);
  void push_ray(int, const T_ray&, const T_cell_tracker&);
  int pop_ray();

  /** This function consolidates all the cache arrays to point to a
      contiguous piece of memory. */
  void consolidate_arrays(array_1&); 
  /** For scalars (monochromatic runs), consolidating is a no-op,
      because the cache variables are scalars and have no dynamic
      storage. */
  void consolidate_arrays(T_float&) {};
  /** For TinyVectors (aux runs), consolidating is a no-op, because
      the cache variablesn have no dynamic storage. */
  template <typename T, int N>
  void consolidate_arrays(blitz::TinyVector<T,N>&) {};

public:
  xfer (/// The grid in which ray tracing is to be done
	T_grid&,
	/// The emission object.
	T_emission&,
	/// The emergence object.
	T_emergence&,
	/// The dust_model object.
	const T_dust_model&,
	/// Seed for the random number generator.
	T_rng_policy& rng,
	/** The lowest intensity for which rays are traced.  Rays of
            lower intensity are dropped. */
	T_float i_ls=1e-4,
	/// The lowest ray intensity for which scatterings are forced.
	T_float i_f=1e-2,
	/// The maximum number of scatterings that are forced.
	T_float n_f= 2,
	/// The maximum intensity before rays are split.
	T_float i_s=1e2,
	/** If true, the xfer object makes a local copy of the
            emergence object and uses this to record emerging
            rays. When the xfer object destructs, the data in the
            local emergence object are added to the object supplied in
            the constructor call.  Using a local copy is advantageous
            since it then becomes unnecessary to use mutexes when
            writing to the emergence object, but of course uses more
            memory. NOTE that currently there are issues with this
            that may cause it not to work, suggest don't use it for
            now. */
	bool keep_local_emergence=false,
	/** If true, radiative intensities are added to the grid. */
	bool add_intensities=true,
	/// The minimum number of scatterings to register.
	int minscat=0,
	/// The maximum number of scatterings to register.
	int maxscat=blitz::huge(int()));
  xfer (const xfer& rhs);
  ~xfer ();
  
  T_rng_policy& rng () {return rng_;};

  // Transfers one ray until it exits the volume or is absorbed
  void shoot(const T_biaser&);
  void shoot_isotropic(const T_biaser&);

  /** Integrates column densities along a ray starting from the
      specified position, going in the specified direction until it
      leaves the grid. This is useful to estimate column densities
      through the grid. */
  T_densities integrate_ray(const vec3d& pos, const vec3d& dir);

  bool uses_local_emergence() const { 
    return local_emergence_.get() ? true : false; };
};

/** Constructor. Sets up the transfer class to handle a specified grid
    and send the output to a specified emergence object. */
template <typename dust_model_type, typename grid_type, typename T_rng_policy>
mcrx::xfer<dust_model_type, grid_type, T_rng_policy>::
xfer(T_grid& grid, T_emission& emi, T_emergence& eme, 
     const T_dust_model& dm, T_rng_policy& rng,
     T_float i_ls, T_float i_f, T_float n_f, T_float i_s,
     bool keep_local_emergence, bool add_intensities,
     int minscat, int maxscat) :
  g (grid), current_cell (), I_last_scatter (i_ls), 
  I_forced (i_f), N_forced (n_f), ray_max_i_(i_s), 
  emergence_ (eme), dust_model(dm), rng_(rng),
  local_emergence_ (keep_local_emergence? new T_emergence (eme): 0),
  emission_ (emi.clone()), 
  dn_(make_thread_local_copy(dust_model.zero_density())),
  add_intensities_(add_intensities),
  n_scat_min_(minscat), n_scat_max_(maxscat)
  // we do not need to initialize the ray here because it has a fully
  // compliant assignment operator that resizes everything as
  // necessary
{
  if (local_emergence_.get()) {
    current_emergence_ = local_emergence_.get();
    local_emergence_->use_lock(false );
  }
  else
    current_emergence_ = &emergence_;

  consolidate_arrays(i_in_);
}


/**  Copy constructor.  Makes a copy of an xfer object, including the
  current cell and random number generator state. */
template <typename dust_model_type, typename grid_type, typename T_rng_policy>
mcrx::xfer<dust_model_type, grid_type, T_rng_policy>::
xfer(const xfer& rhs):
  g (rhs.g), current_cell (rhs.current_cell), 
  I_last_scatter (rhs.I_last_scatter), rng_(rhs.rng_),
  ray_ (rhs.ray_), I_forced (rhs.I_forced), N_forced (rhs.N_forced),
  ray_max_i_(rhs.ray_max_i_), 
  emergence_ (rhs.emergence_), dust_model (rhs.dust_model),
  // if it has a local emergence, we copy that as well
  local_emergence_ (rhs.local_emergence_.get()? 
		   // we shoul copy the LOCAL emergence, right?
		   new T_emergence (*rhs.local_emergence_): 0),
  emission_ (rhs.emission_->clone()), 
  dn_(make_thread_local_copy(dust_model.zero_density())),
  add_intensities_(rhs.add_intensities_),
  n_scat_min_(rhs.n_scat_min_), n_scat_max_(rhs.n_scat_max_)

{
  if (local_emergence_.get()) {
    current_emergence_ = local_emergence_.get();
    local_emergence_->use_lock(false );
  }
  else
    current_emergence_ = &emergence_;

  consolidate_arrays(i_in_);
}

/** Destructor adds the data in the local emergence to the main
    emergence, if applicable. */
template <typename dust_model_type, typename grid_type, typename T_rng_policy>
mcrx::xfer<dust_model_type, grid_type, T_rng_policy>::~xfer ()
{
  // if we've been using a local emergence copy, add it to the
  // global one
  if (local_emergence_.get())
    emergence_ += *local_emergence_; 
};

template <typename dust_model_type, typename grid_type, typename T_rng_policy>
void mcrx::xfer<dust_model_type, grid_type, T_rng_policy>::
consolidate_arrays(array_1&)
{
  const int narr=10;
  // Pad the allocated memory so different threads don't share cache
  // lines.
  const int buffer=2; 
  // allocate the space and set this block to not use locking
  array_2 temp(narr+buffer*2, emission_->zero_lambda().size());
  threadLocal_warn(temp);
  int i=buffer;
  // Because temp is local, we can't return *weak* references.
  norm_.reference(temp(i++,blitz::Range::all()));
  i_in_.reference(temp(i++,blitz::Range::all()));
  i_out_.reference(temp(i++,blitz::Range::all()));
  dtau_.reference(temp(i++,blitz::Range::all()));
  tau_exit_.reference(temp(i++,blitz::Range::all()));
  onemexptau_exit_.reference(temp(i++,blitz::Range::all()));
  scattering_tau_.reference(temp(i++,blitz::Range::all()));
  original_intensity_.reference(temp(i++,blitz::Range::all()));
  intensity_to_camera_.reference(temp(i++,blitz::Range::all()));
  assert(i==narr+buffer-1);
}


/** Generates a ray and propagates it until it's done. This is the
    top-level routine called by the user. Unfortunately it's very
    long. */
template <typename dust_model_type, typename grid_type, typename T_rng_policy>
void
mcrx::xfer<dust_model_type, grid_type, T_rng_policy>::shoot (const T_biaser& bb)
{
  using namespace blitz;

  // copy supplied biaser to the member 
  b = bb;

  hpm::enter_stage(1);
  // Generate a ray
  DEBUG(2,cout << "\nShooting new ray" << endl;);
  emit();

  while(current_cell) {

    DEBUG(2,cout << "ray max I/max L/max bias: " << max(ray_.intensity()) << '\t' << max(ray_.intensity()*norm_) << '\t' << max(b.bias_factor(ray_.intensity())) << endl;);
    DEBUG(1,if(!condition_all(ray_.intensity()< ray_max_i_)) \
	      cout << "WARNING: Max ray intensity " << max(ray_.intensity())<< " at lambda " << maxIndex(ray_.intensity()) << " reflambda " << b.reference_index() << " refintensity " << b.reference1(ray_.intensity()) << endl;);

    // Is ray OK for force-scatter?
    if( (b.reference1(ray_.intensity ())>I_forced) && 
	(ray_.scatterings() < N_forced) ) {
      // Yes - this ray will be force-scattered
      hpm::enter_stage(5);

      // Initialize forced scattering variables. Note that the ray can
      // already have nonzero l at this point if it originated outside
      // the grid. We must take into account that length does NOT
      // refer to length from origin when we use it.
      const vec3d ray_origin = ray_.position();
      traversed.clear();
      traversed.push_back(T_propagation_history (dust_model.zero_density(),
						 ray_.length(), T_cell_tracker()) );
      
      // 1.  Propagate ray to grid boundary without scattering,
      // accumulating optical depth, saving traversed and not adding
      // intensity directly
      hpm::enter_stage(3);
      while(current_cell)
        propagate<false, true, false>();
      hpm::enter_stage(5);

      // 1.5 Add intensities to cells. We can do this before drawing a
      // scattering point because we know that the intensities will
      // decline as exp(-tau) and we'll probably get less noise by
      // doing this analytically.

      // use two temporary arrays for storage and swap them back and forth
      i_in_ = ray_.intensity();

      hpm::enter_stage(4);
      if(add_intensities_)
	for (typename std::vector<T_propagation_history >::const_iterator 
	       i = traversed.begin(), j = i++; i != traversed.end(); ++i, ++j) {
	  assert (i != traversed.begin());
	  assert (i- j == 1);
	  const T_float dl = i->length() - j->length();
	  // optical depth through this cell

	  /// \todo this is actually quite expensive. Can we speed it up?
	  /// Put the difference in dn_ and then shortcut if it's only
	  /// one element?
	  dtau_ = dust_model.total_abs_length( i->n()(tensor::j) - 
					       j->n()(tensor::j) );

	  assert(i_out_.data()!=i_in_.data());
	  i_out_ = i_in_*exp(-dtau_);

	  // The appropriate intensity to add is the length-averaged
	  // mean intensity through the grid cell, which is the
	  // expression below (see notes 080909). As dtau->0, we will
	  // suffer cancellation and possibly division by zero, so we
	  // use expansion for very small dtau.
	  add_intensity (i->cell(), dl, where(dtau_<1e-6,
					      0.5*(i_in_+i_out_),
					      (i_in_-i_out_)/dtau_));
	  swap(i_in_, i_out_);

	}
      hpm::enter_stage(5);
        
      // 2.  Draw an interaction optical depth within the distribution
      // If the medium is empty, there is a degeneracy. NO part of the
      // ray will scatter and we need to deal with this separately below.
      tau_exit_ =  
	dust_model.total_abs_length(ray_.traversed_column()(blitz::tensor::j));
      onemexptau_exit_ = -expm1 (-tau_exit_);
      const T_float onemexpreftau_exit = b.reference1(onemexptau_exit_);

      // we must take care of very small tau_exit values separately to
      // avoid getting poor precision when doing the log(1+x). We use an
      // approximation for very small x, as 0 is not acceptable. (see
      // bug 33).
      const T_float x=-rng_.rnd()*onemexpreftau_exit;
      scattering_reftau_ = x<-1e-5 ? 
	-log(1+x) : -x*(6.+x)/(6.+4.*x);
	
      // 3.1 We need to find out the full scattering_tau vector at all
      // wavelengths.  This is necessary to do before we start adding
      // intensities to the cells, because the ray must be
      // biased. Search the propagation history to find the scattering
      // cell. The iterator will point to the element containing the
      // cell which will scatter, and the n and l at the END of that
      // cell.
      typename std::vector<T_propagation_history >::const_iterator 
	scattering_i = 
	upper_bound(traversed.begin(), traversed.end(), 
		    scattering_reftau_,
		    boost::lambda::_1 < 
		    boost::lambda::bind(&T_dust_model::total_reference_abs_length,
					boost::cref(dust_model),
					boost::lambda::bind(&T_propagation_history::n, 
							    boost::lambda::_2),
					boost::cref(b)));
      assert((scattering_reftau_ > 0) || (scattering_i==traversed.end()));
      const T_cell_tracker& scattering_c(scattering_i->cell());

      // scattering_i will be traversed.end() if the medium is empty!
      // (but not if the medium is *practically* empty)
      if (scattering_i == traversed.end()) {
	// the medium is empty, end the ray
	DEBUG(2, cout << "Forced scattering in empty medium, no scattering." << endl;);
	current_cell.reset();	
	break;
      }

      const T_float endreftau = 
	dust_model.total_reference_abs_length(scattering_i->n(),b);
      const T_float begreftau = 
	dust_model.total_reference_abs_length((scattering_i-1)->n(),b);
      assert(begreftau <= scattering_reftau_);
      assert(endreftau >  scattering_reftau_);

      // Fractional optical depth traversed at the scattering point,
      // which is the same as the fractional column densities
      // traversed
      const T_float fn = 
	(scattering_reftau_ - begreftau)/(endreftau - begreftau);

      // The traversed column densities of the ray at the scattering point
      const T_densities scattering_n 
	( (scattering_i-1)->n() + 
	  (scattering_i->n()-(scattering_i-1)->n()) * fn );

      // Find fractional distance through the cell at the scattering
      // point 
      const T_float dl = scattering_i->length() - (scattering_i-1)->length();
      // (remember that to get to the scattering point we need to
      // subtract the length the ray had gone when it hit the grid)
      const T_float scell_l = 
	(scattering_i-1)->length()-traversed.begin()->length();
      const T_float f = 
	scattering_c.column_to_distance_fraction
	(fn, ray_origin + scell_l*ray_.direction(), ray_.direction(), dl);

      // The path length of the ray from entering the grid to scattering point
      const T_float scattering_l = scell_l + f*dl;

      // The full scattering optical depth vector
      scattering_tau_ =
	dust_model.total_abs_length(scattering_n(blitz::tensor::j));

      // 3.2 Calculate bias factor due to optical depth of scattering
      // based on the forced probabilities. 
      DEBUG(1,T_lambda bias(b.bias_factor(scattering_tau_*exp(-scattering_tau_)/(onemexptau_exit_))); \
	    if(!condition_all(bias<ray_max_i_)) \
	      cout << "WARNING: Large bias factor  " << max(bias)<< " at lambda " << maxIndex(bias) << " reflambda " << b.reference_index() << " scattering reftau " << scattering_reftau_ << " reftau_exit " << b.reference1(tau_exit_) << " at "<<__FILE__ <<':'<<__LINE__ <<endl;);

      // we need to deal with very small scattering_tau separately
      // because we get cancellation and 0/0. For tau_exit<1e-6, we
      // use the taylor expansion of the expression, which is
      // scattering_tau/tau_exit. This gives a max approximation error
      // of about 1e-6, which is way below our other errors.
      ray_.bias(b.bias_factor
		(where(tau_exit_<1e-6, 
		       scattering_tau_/tau_exit_, 
		       scattering_tau_*exp(-scattering_tau_)/(onemexptau_exit_))));

      // save intensity so we can restore it below
      original_intensity_ = ray_.intensity();
	
      // Check for degenerate case where the ENTIRE ray exited and it
      // doesn't make sense to deal with the forced scattering ray and
      // we bail out and the ray is done.  this was already checked
      // for above so should never happe
      assert(scattering_l>0);

      // 4. update the ray to the interaction point, and scatter it
      ray_.set_position(ray_origin);
      // set ray intensity to part that will scatter
      ray_.bias(onemexptau_exit_);
      ray_.propagate(scattering_l, scattering_n);
      current_cell = scattering_c;
      assert (current_cell);

      // This assertion is not quite safe because ray_origin can be on
      // the grid boundary in which case locate can give zero.  In
      // cases of suddenly very high optical depth in cells, where the
      // ray will scatter very close to the cell face, numerical
      // effects can also make the ray be in the immediately preceding
      // cell. This is more likely to happen if the emission point is
      // far away, as the cells will be small dl's on top of a large
      // value just to make it to the grid. can we reformulate the
      // math here to be less sensitive to this?
      current_cell.assert_contains(ray_.position());

      scatter ();

      // This completes this ray segment, we now start over and see
      // whether we should force scatter again

    } // if force scatter
    else {
      hpm::enter_stage(0);

      // The ray is not forced to scatter, we just need to draw
      // scattering tau from the un-forced distribution
      scattering_reftau_ = - log (rng_.rnd ());

      traversed.clear();
      traversed.push_back(T_propagation_history (dust_model.zero_density(),
						 ray_.length(), T_cell_tracker()) );
      
      // Propagate ray until scatter or exit, save traversed
      // intensities, but do not add to cells directly
      hpm::enter_stage(3);
      while( (!propagate<true, true, false>()) && current_cell);
      hpm::enter_stage(0);

      // calculate bias factor due to optical depth of scattering
      // based on the unforced probabilities. We need to bias BEFORE
      // adding the intensities because the intensities to add depend
      // on whether the ray exited or not.
      scattering_tau_ =
	dust_model.total_abs_length(ray_.traversed_column()(blitz::tensor::j));

      // the bias depends on whether the ray exited the cell or not
      if (current_cell) {
	// bias factor for rays that did not exit is exp(-tau)/tau.
	DEBUG(1,T_lambda bias(b.bias_factor(exp(-scattering_tau_)* \
					    scattering_tau_));	  \
	      if(!condition_all(bias<ray_max_i_))		  \
		cout << "WARNING: Large bias factor  " << max(bias)<< " at lambda " << maxIndex(bias) << " reflambda " << b.reference_index() << " scattering reftau " << scattering_reftau_ << " at " << __FILE__ <<':' <<__LINE__ <<endl;);
	ray_.bias(b.bias_factor(exp(-scattering_tau_)*scattering_tau_));
      }
      else {
	// bias factor for rays that exit is exp(-tau)
	DEBUG(1,T_lambda bias(b.bias_factor(exp(-scattering_tau_))); \
	      if(!condition_all(bias<ray_max_i_))		    \
		cout << "WARNING: Large bias factor  " << max(bias)<< " at lambda " << maxIndex(bias) << " reflambda " << b.reference_index() << " scattering reftau " << scattering_reftau_ << " at " << __FILE__ <<':' <<__LINE__ <<endl;);

	ray_.bias(b.bias_factor(exp(-scattering_tau_)));
      }

      // add intensities now that we have biased the ray
      hpm::enter_stage(4);
      if(add_intensities_)
	for (typename std::vector<T_propagation_history >::const_iterator 
	       i = traversed.begin(), j = i++; i != traversed.end(); ++i, ++j) {
	  assert (i != traversed.begin());
	  assert (i- j == 1);
	  const T_float dl = i->length() - j->length();
	  add_intensity (i->cell(), dl);
	}
      hpm::enter_stage(0);

      if (current_cell) {
	// if current_cell is set, we exited the loop above because
	// propagate returned true indicating that the ray should be
	// scattered (we do it here to avoid recursive calling of
	// propagate)
	scatter();
      }
      else {
	DEBUG(2,cout << "Ray left volume" << endl;);
      }
    } // if not force scatter

    // coming out here, we have propagated the ray to either the
    // scattering point or until it left the grid

    // if the ray didn't leave, check if it should be split
    if(current_cell)
      check_ray_split();
    
    // If the ray lost russian roulette or the ray left, current_cell
    // is set to zero by scatter.  If so, the ray is done. 

    // We must now check the stack to see if there are split rays to
    // process. Popping them resets current cell, so in that case we
    // continue looping.
    if(!current_cell && !ray_stack_.empty()) {
      const int r = pop_ray ();
      if(r >= 0)
	// true means the ray was split pre-scatter and must be scattered now
	scatter(r);
    }
  } // while current_cell

  // Coming out here, the ray is either below our intensity cutoff, or
  // it left the grid.  In any case, we are done with this ray.
}


/** Pops a ray, including its current cell and pre/post scatter
    status, off of the ray stack and makes it the current ray.  */
template <typename dust_model_type, typename grid_type, typename T_rng_policy>
void
mcrx::xfer<dust_model_type, grid_type, T_rng_policy>::
push_ray(int rw,
	 const T_ray& r, 
	 const T_cell_tracker& c) 
{
  DEBUG (2, cout << "Pushing ray with max intensity " << max(r.intensity()) << endl;);
  ray_stack_.push_back(ray_stack_object<T_ray, T_cell_tracker>(rw, r, c));
}


/** Pops a ray, including its current cell and pre/post scatter
    status, off of the ray stack.  */
template <typename dust_model_type, typename grid_type, typename T_rng_policy>
int mcrx::xfer<dust_model_type, grid_type, T_rng_policy>::pop_ray()
{
  ray_ = ray_stack_.back().ray_;
  current_cell = ray_stack_.back().cell_;
  const int r = ray_stack_.back().resume_with_;
  DEBUG (2, cout << "Popping " << (r<0?"post":"pre") << "scatter ray with max intensity " << max(ray_.intensity()) << endl;);
  
  ray_stack_.pop_back();
  return r;
}

/** Check ray max intensity and possibly split. */
template <typename dust_model_type, typename grid_type, typename T_rng_policy>
void
mcrx::xfer<dust_model_type, grid_type, T_rng_policy>::check_ray_split(int r)
{
  ASSERT_ALL(ray_.intensity() == ray_.intensity());

  if ((ray_max_i_ == 0) || !current_cell) 
    // no max i set, check unnecessary
    return;

  const T_float maxi = max(ray_.intensity());
  DEBUG(3,cout << "Checking ray intensity: max I " << maxi << endl;);
  if(maxi > ray_max_i_) {
    // ray should be split
    const int n = int(ceil(maxi/ray_max_i_));
    DEBUG(1,cout << "Splitting " << (r<0?"post":"pre") << "scatter ray with max intensity " << maxi << " into " << n << endl;);

    ray_.bias(1./n);
    for(int i=0;i<n-1;++i)
      // we push n-1 onto the stack, because the current ray is still current
      push_ray(r, ray_, current_cell);
  }
}


/** Emits a ray according to the emission distribution, and calculates
    the direct flux to the cameras. */
template <typename dust_model_type, typename grid_type, typename T_rng_policy>
void mcrx::xfer<dust_model_type, grid_type, T_rng_policy>::emit()
{
  // We call on the emission object to emit the ray and give us the
  // normalization, velocity, and angular distribution of emission.
  ray_distribution_pair<T_lambda> e = emission_->emit(b, rng_);
  ray_ = *e.ray_pointer();
  norm_ = e.normalization();

  /// \todo this needs to be evaluated.
  /// if the normalization is zero, we can safely set the intensity to
  /// be zero, too. this eliminates any chance of getting high bias
  /// factors at those wavelengths
  //ray_.set_intensity(where(norm_>0.0, ray_.intensity(), 0.0));

  ray_.set_traversed_column(dust_model.zero_density() );

  assert (condition_all(ray_.intensity() >= 0));    

  // Emerge the direct radiation getting to the cameras
  emerge_from_point (e.distribution(), e.velocity());

  // now set the doppler shift of the ray based on particle velocity and ray
  // direction (can't do this before the direct radiation has been treated.)
  const T_float d = 1+dot(e.velocity(), ray_.direction())/constants::c0;
  ray_.add_doppler(d);
  DEBUG(1,std::cout << "Dopplerfactor " << d << " in emit"<< std::endl;);

  // propagate the ray to the grid boundary, if it is outside, and
  // locate the current cell
  propagate_external ();
}


/** Propagates the ray from its position to the end of the current
    grid cell or until an interaction takes place.  If do_scatter is
    false, we ignore scattering and just keep on propagating until we
    hit the boundary or max_length. (This is used for forced
    scattering.) If keep_traversed is true, the traversed cells are
    added to the propagation_history vector, if it's false the
    intensities are added to the cells here. Because these two flags
    are only used in a static sense, they are made template parameters
    so the compiler can do the appropriate optimization. */
template <typename dust_model_type, typename grid_type, typename T_rng_policy>
template <bool do_scatter, bool keep_traversed, bool add_intensity_here>
bool
mcrx::xfer<dust_model_type, grid_type, T_rng_policy>::
propagate (T_float max_length)
{
  assert (current_cell);
  assert (max_length >= ray_.length());

  // save the current cell so we can access it.
  const T_cell_tracker cc(current_cell);

  // Find the cell boundary in the direction of travel. NOTE THAT for
  // efficiency's sake is THIS UPDATES current_cell and dn_. If you
  // want the cell that we entered in, use the cc tracker object.

  // (this routine also makes sure we don't propagate past max_length,
  // resets current_cell if we hit max_length, and ensures the column
  // density returned in dn_ is appropriate in case we stop due to max
  // length)

  const T_float dl
    (current_cell.intersection_from_within(ray_, max_length, dn_));
  
  assert (dl > -1e-10);

  // Check if we have reached optical depth for absorption/scattering
  const T_float reftau = 
    dust_model.total_reference_abs_length(ray_.traversed_column (), b);
  const T_float dreftau = dust_model.total_reference_abs_length(dn_, b);
  assert (dreftau >= -1e-10);

  // true if scattering should happen in this grid cell
  const bool scatter_now = 
    do_scatter && (scattering_reftau_ < reftau + dreftau);
  
  // 2 things can happen: scattering or just moving to the edge of the
  // grid cell.  
  if(scatter_now) {
    // Scattering is enabled and is to take place in this grid cell

    // column density fraction to scattering event
    const T_float fn = (scattering_reftau_ - reftau)/dreftau;

    // Calculate fractional distance to absorption/scattering event
    // based on fractional column density (this needs to be done by
    // the cell_tracker since the cells may not have uniform density).
    const T_float fl = 
      cc.column_to_distance_fraction(fn ,ray_.position(), ray_.direction(), dl);

    assert (fl < 1);
    // This assertion can fail if the density in the cell is so high that
    // the scattering distance underflows to 0. This usually means the density
    // in the cell is catastrophically wrong, like wrong units.
    assert (fl > 0);

    // Update ray position to pos of scattering event
    ray_.propagate (fl*dl, fn*dn_);
  }
  else {
    // Either scattering is not turned on or it's not time for
    // scattering. In both cases we just propagate to the edge of the
    // cell.

    ray_.propagate(dl, dn_);      
  }

  if (keep_traversed)
    // Update the propagation history
    traversed.push_back(T_propagation_history 
			(independent_copy(ray_.traversed_column()),
			 ray_.length(),
			 cc));
  if (add_intensity_here)
    // add intensity to the cell
    add_intensity (cc, dl);

  if (scatter_now) {
    // we scattered in this cell, which means we are still in the cell
    // we entered. In this case we need to roll back current_cell,
    // which was automatically set to the next cell by
    // intersection_from_within
    current_cell = cc;
  }
  
  return scatter_now; 
}


/** Scatters the ray.  The ray is scattered at the current position by
    calling the scatter function of the absorber in the current cell.
    After this, flux emerging to the cameras from the scattering is
    calculated.  If the ray is completely absorbed (too low intensity)
    after the scattering, current_cell is set to 0 to indicate that
    this ray is done. */
template <typename dust_model_type, typename grid_type, typename T_rng_policy>
void
mcrx::xfer<dust_model_type, grid_type, T_rng_policy>::scatter(int resume_with)
{
  assert (current_cell);
  DEBUG(2,cout << "Scattering ray for " << ray_.scatterings() << " time" << endl;);

  int i;
  if (resume_with < 0) {
    // Get scatterer that will do the interaction from dust_model. This
    // includes biasing due to different opacities.
    i = 
      dust_model.draw_scatterer(current_cell.data()->get_absorber().densities(),
				b,
				ray_,
				rng_);
    // attenuate ray intensity based on albedo. This is where
    // luminosity gets absorbed into the grid. In debugging mode, we
    // keep track of this for energy conservation checks.
#ifdef MCRX_DEBUG_LEVEL
  const T_lambda pre_intensity(independent_copy(ray_.intensity()));
#endif

    assert (condition_all(ray_.intensity() >= 0));    
    ray_.bias(dust_model.get_scatterer(i).albedo());

    // we only do this when debugging because it requires mutex
    // locking and isn't used except to test for energy conservation.
    DEBUG(1,dust_model.add_absorbed					\
	  (T_lambda((pre_intensity - ray_.intensity())*norm_)););
    DEBUG(1,current_cell.data()->add_absorbed		\
	  ((pre_intensity - ray_.intensity())*norm_););
    
    // At this point the ray has the intensity that is outgoing from
    // the grain, but yet unaffected by outgoing direction. This is
    // the appropriate intensity to be sent to cameras. Note that we
    // have NOT drawn a scattering direction yet, and
    // emerge_from_point restores the ray direction. Conceptually, the
    // ray has already scattered, so we manually increase the nscat
    // counter so it represents correctly the light after scattering.
    ++ray_.scatterings();
    emerge_from_point (dust_model.get_scatterer(i).phase_function(),
		       current_cell.data()->get_absorber().velocity());
    // and set it back
    --ray_.scatterings();
    
    check_ray_split(i);
  }
  else {
    // resume with was set, so this ray was split after emerging.
    // Then we resume with the scattering part.
    i = resume_with;
  }
  
  {
    // save the incoming direction for doppler calculation
    const vec3d indir(ray_.direction());
    
    // Now we ask the grain to scatter the ray in a direction and apply
    // biasing due to outgoing direction. We don't care about the
    // returned distribution now.
    dust_model.get_scatterer(i).scatter(b, ray_, rng_);
    
    // Apply doppler shift for scattering event
    const T_float d = 1.0 + dot(current_cell.data()->get_absorber().velocity(),
				ray_.direction() - indir)/constants::c0;
    ray_.add_doppler(d);
    DEBUG(1,cout << "Dopplerfactor " << d << ", velocity " << current_cell.data()->get_absorber().velocity() << " in scatter"<< endl;);
  }
  current_cell.reset_intersections();

  // This completes the scattering event.

  // Do Russian Roulette for low-intensity rays

  // Test if the max intensity is below the RR threshold
  if (max(ray_.intensity()) < I_last_scatter) {
    
    // If ray has zero intensity it doesn't matter if it passed, it's dead
    if(condition_all(ray_.intensity() == 0)) {
      DEBUG(2,cout << "Ray intensity ==0, killing"<<endl;);
      current_cell.reset();
    }
    
    DEBUG(2,cout << "Ray max intensity " << max(ray_.intensity()) << ",  playing Russian Roulette... ";);
    const double P_rr=0.1; 
    if(rng_.rnd() < P_rr) {
      // survived, boost intensity by probability factor
      DEBUG(2,cout <<"Won" << endl;);
      ray_.bias(1./P_rr);
    }
    else {
      // lost, kill it. This is done by resetting current_cell.
      DEBUG(2,cout <<"Lost" << endl;);
      current_cell.reset();
    }
  }
}


/** Calculates the radiation that would reach the cameras from the
    current ray position as a result of emission or scattering. */
template <typename dust_model_type, typename grid_type, typename T_rng_policy>
void
mcrx::xfer<dust_model_type, grid_type, T_rng_policy>::
emerge_from_point (const angular_distribution<T_lambda>& d,
		   const vec3d& velocity)
{
  assert (condition_all(ray_.intensity() >= 0));   
  const int curstage=hpm::current_stage();
  hpm::enter_stage(6);

  // Save these so we can restore them when we are done
  const T_cell_tracker cc = current_cell;
  const vec3d original_direction = ray_.direction();
  const T_float original_doppler = ray_.doppler();
  const vec3d ray_origin = ray_.position();

  // go through the cameras one by one
  typename T_emergence::iterator stop= current_emergence_->end();
  for (typename T_emergence::iterator c = current_emergence_->begin();
       c != stop; ++c) {

    // 1. Find the direction to the camera c 
    ray_.set_position(ray_origin );
    ray_.set_doppler(original_doppler);
    const vec3d delta = (*c)->position () - ray_.position();
    const T_float dist2 = dot (delta, delta);
    const T_float dist = sqrt(dist2);
    const T_float idist2 = 1./dist2;
    ray_.set_direction(delta/dist );
    current_cell.reset_intersections();

    // 1.1 assign doppler
    const T_float doppler = 1.0 + dot(velocity, 
				      ray_.direction() - original_direction)/
      constants::c0;
    ray_.add_doppler(doppler);
    DEBUG(1,cout << "Dopplerfactor " << doppler << " in emerge from point"<< endl;);

    // 2.  See if the ray is outside the grid volume, and in that
    // case, if it will hit the grid volume (not going past the camera)
    if (!current_cell)
      propagate_external (dist);

    // 3.  Assuming it hit the grid, propagate the ray through the
    // grid (up to max length=dist), without scattering or
    // saving/adding intensities.
    hpm::enter_stage(3);
    while (current_cell)
      propagate<false, false, false> (dist);
    hpm::enter_stage(7);

    // 4.  Now add the flux that would get to the camera to the camera.
    
    // first write the distribution function in place into the itc
    // array to avoid allocating
    d.distribution_function(ray_.direction(), 
			    original_direction,
			    intensity_to_camera_);
    // calculate the optical depth and put it in dtau_ (we want to put
    // this result in a temporary because it is an index expression
    // and hence not vectorized.  Otherwise I also think we redo the
    // sum for each lambda element)
    hpm::enter_stage(8);
    dtau_ = dust_model.total_abs_length(ray_.traversed_column()(blitz::tensor::j)); 
    // and now multiply everything together
    hpm::enter_stage(9);
    intensity_to_camera_ *= ray_.intensity() * norm_ * idist2 * exp (-dtau_);
    //intensity_to_camera_ *= ray_.intensity() * norm_ * idist2 * exp (-dust_model.total_abs_length(ray_.traversed_column()(blitz::tensor::j)));
   
    assert (condition_all(intensity_to_camera_ >= 0));    
    DEBUG(2,cout << "Adding ray to camera, max i " \
	  << max(ray_.intensity()* \
		 d.distribution_function(ray_.direction(), original_direction))\
	  << endl;);
    DEBUG(1, array_1 temp(ray_.intensity()* \
			  d.distribution_function(ray_.direction(), \
						  original_direction)); \
	  if( max(temp) > ray_max_i_ ){ \
	    cout << "WARNING: Adding ray to camera, max intensity " \
		 << max(temp) << " at lambda " << maxIndex(temp) \
		 << " reflambda " << b.reference_index() \
		 << endl;};);
    hpm::enter_stage(2);
    add_to_camera(*c, ray_origin, intensity_to_camera_, ray_.doppler());
    hpm::enter_stage(6);
    
    current_cell = cc; 
  }

  // Restore ray to the incoming state
  ray_.set_position(ray_origin );
  ray_.set_direction(original_direction);
  ray_.set_doppler(original_doppler);
  current_cell.reset_intersections();

  hpm::enter_stage(curstage);
}

      
/** Propagates a ray outside the grid to the grid boundary, if it hits
    the grid.  This is used if the ray is emitted outside the grid. */
template <typename dust_model_type, typename grid_type, typename T_rng_policy>
void 
mcrx::xfer<dust_model_type, grid_type, T_rng_policy>::
propagate_external (T_float max_length)
{
  assert (!current_cell);

  // ask the grid whether the ray will enter a valid cell
  const std::pair<T_cell_tracker, T_float> tl = 
    g.intersection_from_without(ray_ );

  current_cell = tl.first;

  if (current_cell && tl.second < max_length) {
    // the ray will enter a valid cell within max_length

    ray_.propagate(tl.second, dust_model.zero_density());
  }

  assert(ray_.length()<=max_length);
}


/** Emits a ray and immediately adds the radiation to the cameras
    without scattering, as if in a vacuum. This is the "budget
    version" of shoot and is used when we just want to make an image
    of the emission without doing any radiative transfer. */
template <typename dust_model_type, typename grid_type, typename T_rng_policy>
void mcrx::xfer<dust_model_type, grid_type, T_rng_policy>::shoot_isotropic (const T_biaser& bb)
{
  // copy supplied biaser to the member 
  b = bb;

  // We call on the emission object to emit the ray (the peculiar
  // maneuver is needed because we need to return the reference to the
  // angular_distribution object as well as the emitted ray)
  ray_distribution_pair<T_lambda> e = emission_->emit(b, rng_);
  ray_ = *e.ray_pointer();
  norm_ = e.normalization();
  const T_float original_doppler = ray_.doppler();

  // go through the cameras one by one
  typename T_emergence::iterator stop= current_emergence_->end();
  for (typename T_emergence::iterator c = current_emergence_->begin();
       c != stop; ++c) {

    // 1. Find the direction to the camera c 
    const vec3d delta = (*c)->position () - ray_.position();
    const T_float dist2 = dot (delta, delta);
    const T_float dist = sqrt(dist2);
    ray_.set_direction(delta/dist );

    // 1.1 assign doppler
    ray_.set_doppler(original_doppler);
    const T_float doppler = 1.0 + dot(e.velocity(), ray_.direction())/
      constants::c0;
    ray_.add_doppler(doppler);
    DEBUG(1,cout << "Dopplerfactor " << doppler << " in shoot_isotropic"<< endl;);

    // 2.  Add the flux to the camera
    const T_lambda intensity_to_camera 
      (make_thread_local_copy
       (ray_.intensity()*norm_*
	e.distribution().distribution_function(ray_.direction(), vec3d (0,0,0))
	/ dist2 ));
    DEBUG(3,std::cout << "Adding Ratio " << intensity_to_camera(535)/intensity_to_camera(532) << ' '<< intensity_to_camera(615)/intensity_to_camera(614) << endl;);

    assert (condition_all(intensity_to_camera >= 0.));    
    (*c)->add(ray_.position(), intensity_to_camera, ray_.doppler());
  }
}


/** Add intensity to the cell. */
template <typename dust_model_type, typename grid_type, typename T_rng_policy>
void 
mcrx::xfer<dust_model_type, grid_type, T_rng_policy>::
add_intensity (const T_cell_tracker& c,
	       T_float dl)
{
  add_intensity(c, dl, ray_.intensity());
}


/** Add intensity to the cell. */
template <typename dust_model_type, typename grid_type, typename T_rng_policy>
template <typename T>
void 
mcrx::xfer<dust_model_type, grid_type, T_rng_policy>::
add_intensity (const T_cell_tracker& c,
	       T_float dl,
	       const blitz::ETBase<T>& intensity)
{
  if(!add_intensities_) 
    return;

  const T& intens=intensity.unwrap();

  assert(c);
  ASSERT_ALL(dl*norm_*intens<1e60);

  DEBUG(3,cout << "Adding intensity max: " << max(dl*norm_*intens) << endl;);

  c.data()->get_absorber().add_intensity (dl*norm_*intens);
}

/** Integrates the column densities along the ray. */
template <typename dust_model_type, typename grid_type, typename T_rng_policy>
mcrx::T_densities 
mcrx::xfer<dust_model_type, grid_type, T_rng_policy>::
integrate_ray(const vec3d& pos, const vec3d& dir)
{
  ray_=T_ray(pos, dir, dust_model.zero_density(), T_lambda());
  current_cell = g.locate( ray_.position() );
  if (!current_cell) {
    // The ray is outside the grid, we need to see if it
    // will hit and propagate it to that point (it also sets current_cell)
    propagate_external ();
  }
  // propagate it through the grid until it exits
  while(current_cell)
    propagate<false, false, false>();
  
  // return column density
  return ray_.traversed_column();
}


#endif


